#!/usr/bin/env python3

import numpy as np
import numpy.ma as ma
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import cartopy.crs as ccrs
import sys

lon = 0
lat = 0
# proj = ccrs.Orthographic(central_longitude=lon, central_latitude=lat)
proj = ccrs.PlateCarree()
proj._threshold /= 100

plt.figure(figsize=(8, 8))
ax = plt.axes(projection=proj)
ax.coastlines(resolution='110m', lw=0.1)
ax.set_global()

mesh = xr.open_dataset(sys.argv[1])
cell_lon = mesh.cell_lon.to_masked_array()
cell_lat = mesh.cell_lat.to_masked_array()
vtx_lon  = mesh.vtx_lon .to_masked_array()
vtx_lat  = mesh.vtx_lat .to_masked_array()

ax.scatter(cell_lon, cell_lat, s=0.5, c='r', transform=ccrs.PlateCarree())
ax.scatter(vtx_lon , vtx_lat , s=0.5, c='g', transform=ccrs.PlateCarree())

for k in range(6):
	for j in range(mesh.y.size):
		for i in range(mesh.x.size):
			if ma.is_masked(cell_lon[k,j,i]): continue
			lon = [vtx_lon[k,j,i], vtx_lon[k,j,i+1], vtx_lon[k,j+1,i+1], vtx_lon[k,j+1,i], vtx_lon[k,j,i]]
			lat = [vtx_lat[k,j,i], vtx_lat[k,j,i+1], vtx_lat[k,j+1,i+1], vtx_lat[k,j+1,i], vtx_lat[k,j,i]]
			ax.plot(lon, lat, c='b', lw=0.3, transform=ccrs.Geodetic())

# Plot domain outlines.
outlines = np.array([
	[[315, -35.2643896827547], [45 , -35.2643896827547]],
	[[45 , -35.2643896827547], [45 ,  35.2643896827547]],
	[[45 ,  35.2643896827547], [315,  35.2643896827547]],
	[[315,  35.2643896827547], [315, -35.2643896827547]],
	[[45 , -35.2643896827547], [135, -35.2643896827547]],
	[[135, -35.2643896827547], [135,  35.2643896827547]],
	[[135,  35.2643896827547], [45 ,  35.2643896827547]],
	[[135, -35.2643896827547], [225, -35.2643896827547]],
	[[225, -35.2643896827547], [225,  35.2643896827547]],
	[[225,  35.2643896827547], [135,  35.2643896827547]],
	[[225, -35.2643896827547], [315, -35.2643896827547]],
	[[315,  35.2643896827547], [225,  35.2643896827547]]
])
for line in outlines:
 	ax.plot(line[:,0], line[:,1], c='black', lw=0.6, transform=ccrs.Geodetic())

plt.show()
